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ABSTRACT 

The dynamical system studied in previous papers of this series, namely a bound time-like 
geodesic motion in the exact static and axially symmetric space-time of an (originally) 
Schwarzschild black hole surrounded by a thin disc or ring, is considered to test whether the 
often employed “pseudo-Newtonian” approach (resorting to Newtonian dynamics in gravita¬ 
tional potentials modified to mimic the black-hole field) can reproduce phase-space properties 
observed in the relativistic treatment. By plotting Poincare surfaces of section and using two 
recurrence methods for similar situations as in the relativistic case, we find similar tendencies 
in the evolution of the phase portrait with parameters (mainly with mass of the disc/ring and 
with energy of the orbiters), namely those characteristic to weakly non-integrable systems. 
More specifically, this is true for the Paczyriski-Wiita and a newly suggested logarithmic po¬ 
tential, whereas the Nowak-Wagoner potential leads to a different picture. The potentials and 
the exact relativistic system clearly differ in delimitation of the phase-space domain accessi¬ 
ble to a given set of particles, though this mainly affects the chaotic sea whereas not so much 
the occurrence and succession of discrete dynamical features (resonances). In the pseudo- 
Newtonian systems, the particular dynamical features generally occur for slightly smaller 
values of the perturbation parameters than in the relativistic system, so one may say that the 
pseudo-Newtonian systems are slightly more prone to instability. We also add remarks on nu¬ 
merics (a different code is used than in previous papers), on the resemblance of dependence of 
the dynamics on perturbing mass and on orbital energy, on the difference between the New¬ 
tonian and relativistic Bach-Weyl rings, and on the relation between Poincare sections and 
orbital shapes within the meridional plane. 
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1 INTRODUCTION 

Newton’s theory of gravity is still being used in treating many as- 
trophysical systems, because general relativity is (i) often not nec¬ 
essary in weak-field problems, while (ii) often practically inappli¬ 
cable (or only applicable numerically) in strong-field ones. Under 
both circumstances, various approximation methods have been de¬ 
veloped, including, above others, “linearized theory of gravity” and 
post-Newtonian or post-Minkowskian expansions, as well as ad 
hoc effective descriptions like those based on “pseudo-Newtonian” 
potentials. The well-justified small-parameter expansions are typi¬ 
cally reliable in weak-field cases, but in strong field they are inaccu¬ 
rate unless brought to higher expansion orders. The ad hoc formu¬ 
las, though not derived by any sound approximation scheme, may 
be quite simple yet still work well in strong field, but much caution 
is in place, because they often mimic certain features of the prob- 
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lem accurately, while badly misrepresenting the others. Depending 
on particular approach, it may be difficult to specify which kinds of 
errors and of what sizes it brings, the more so if one does not know 
the stability properties of the exact general relativistic solution or if 
such a solution is not even available at all. 

One of thorough ways to assess the practical quality of a given 
description of a given gravitational field, or at least its general dif¬ 
ference from another description, is to study the motion of free 
test particles by methods used in the theory of dynamical systems. 
Though it is problematic to directly compare trajectories of differ¬ 
ent dynamical systems and hence to quantify their relative devia¬ 
tion, it is still possible to compare the systems’ overall “dynami¬ 
cal portraits” and the latter’s dependence on parameters. Needless 
to say, the same methods can reveal the effect of various physical 
perturbations imposed on a given system within the same theory 
or approximation; similarly, they can also be employed to test and 
compare numerical codes. 

_ I n the previous three papers of th is series dSemerak & Suko^ 

|2010ll^ ^ ISukova & Semerakll2013h (below referred to as papers 
I, II and III, respectively), we studied the field of a Schwarzschild- 
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like black hole surrounded by a concentric thin disc or ring, as 
described hy exact static and axially symmetric solution of Ein¬ 
stein equations. Motivated hy astrophysical black holes surrounded 
by accretion (or galactic “circumnuclear”) structures, we analysed 
the gravitational influence of the additional matter on a long-term 
dynamics of time-like geodesics and showed, hy several different 
methods, that it can make the dynamics chaotic. In the present pa¬ 
per we compare the previous relativistic results with a similar anal¬ 
ysis carried out within pseudo-Newtonian description. More specif¬ 
ically, we emulate the Schwarzschild gravitational field by several 
simple “pseudo-Newtonian” potentials, while the disc or ring are 
described by their Newtonian potential (which equals the first of 
the two metric functions appearing in the relativistic description). 

Besides describing the gravitational field and the free test- 
particle motion in a different way, we also use a different numerical 
code to follow the trajectories: whereas the relativistic geodesic- 
equation system was solved, in previous papers, by the Runge- 
Kutta (or rather the Hufa) 6th-order method with variable proper¬ 
time step, here we solve the Newton ian equations of m otion by 
appropriate geometric integrators (see lHairer et al.ll2006h . specifi- 
ca lly in the thin-disc case we have develop ed an integrator inspired 
bv ISevrich & Lukes-GerakopoulosI ( 12012h and endowed with a spe¬ 
cial treatment of the field jump across the disc. In spite of these 
significant differences, we have arrived at a similar dynamical pic¬ 
ture, which justifies the observations made in either of the ways. 
However, there still occur differences with respect to the exact 
Schwarzschild picture, and mainly between the different pseudo- 
Newtonian potentials; some of the latter even do not seem to be 
reasonably applicable. 

After a short note on previous results that have appeared in 
the literature, we specify the pseudo-Newtonian description of our 
gravitational fields in section and review basic properties of mo¬ 
tion in their backgrounds in section Then in section |4] we give 
a basic information about the codes employed. The main section 
[ 5 ] brings the comparison between exact relativistic and pseudo- 
Newtonian results, using Poincare surfaces of section and two re¬ 
currence methods. We add there special notes i) on the link between 
the dependence on perturbing mass and on orbital energy; ii) on 
a different character of the relativistic Bach-Weyl ring and of its 
Newtonian counter-part; and iii) we also point out (and illustrate) 
that the Poincare sections represent only partial information about 
the orbits. Concluding remarks then close the paper. 


1.1 Previous results from the literature 

A similar system we co nsider here was studied by 
IVokrouhlickv & Karas! (Il998h within Newtonian description 
and with motivation stemming from a long-term evolution of 
stars orbiting the black holes (with accretion discs) in galac¬ 
tic nuclei. The authors represented the central body by the 
—Mjr potential and the thin disc by the Kuzmin potential 
—M/+ [A + \z\)'^, denoting by A4 the disc mass, by p 
and 2 ; cylindrical coordinates and by A a free constant, while 
also taking into account mechanical effect of the disc on the test 
orbiter (hydrodynamical drag). The main conclusion was that 
“any consistent model of the star-disc interaction has to take 
the influence of the disc gravity into account, in addition to the 
effects of direct collisions with gaseous material”. The long-term 
dynamics was found to be sensitive to a particular model of the 
disc, especially to the radial profile of its surface density, whereas 
much less to the total mass of the disc. 

The pseudo-'Newtoma.n potentials have been employed in 


many papers on accretion flows around black holes, but only a few 
times in stu dying the chaotic regimes o f motion in perturbed black- 
hole fields. iGueron & Letelien (1200 ih compared the free-motion 
dynamics around a Schwarzschild black hole and around a New¬ 
tonian point centre, when superposed with a dipolar field. They ob¬ 
served that the black-hole system became more chaotic (than the 
exact case) when the centre was simulated by the Paczyhski-Wiita 
pseudo-potentia l, mainly if incorpo rating special relativistic equa¬ 
tion of motion. ISelaru et al.l i2005h studied the Newtonian circu¬ 
lar Hill’s restricted three-body problem while describing the pri- 
mary by the Schwarzschild -type potential Ajr + Bfr^. Similarly, 
ISteklain & Leteli^ ( l2006h compared the Hill problem involving 
the Paczyhski-Wiita pseudo-potential with the original Newtonian 
version, concluding that the pseudo-Newtonian case is usually - 
but not always - more unstable than its Newtonian counter-part. 

Several papers have also tried to incorporate, within the 
pseudo-Ne wtonian approach, draggin g effects due to rotation of 
the centre. ISteklain & Leteli^ ( I2OO9I1 thus found that the orbits 
counter-rotating with respect to the centre are more unstable than 
the co-rotating ones. IWang & Wul ( 1201 ih superposed a rotating 
“pseudo black hole” with a quadrupole halo in order to analyse the 
emission of gravitational waves from orbiting particles; the radi¬ 
ated amplitude and power were observed to be c losely related to the 
degree of orbital chaoticity. The same authors dWang & Wull20l3) 
also used their model in order to discuss how the geodesic dynam¬ 
ics responds to the centre’s spin and to quadrupole perturbation; 
they found, in particular, that the centre’s rotation rather attenuates 
the instability. The dynamics of charged particles in the field of a 
magnetized compact obj ect described in a p seudo-Newtonian man¬ 
ner was then studied by IWang et ^ ( 1201 3h and instabilities were 
identified using the “fast Lyapunov indicator”. 

The advance to the pseudo-Newto nian imitation o f spinn ing 
fields mainly followed the proposal by lArtemova et al.l ( Il996h of 
two simple potentials for the Kerr black hole. Recently, these 
have been checked against a slightly different formula (as well as 
against the “benchmark” of the Paczyhski- Wiita potential) on the 
behav iour of circular-orbit acceleration by iKaras & Abramowiczl 
( l2015h . A m ore elaborate pseudo-N e wtoni a n “fit” of Kerr was 
presen ted by IChakrabarti & Mondall ( l2006l) . Ilvanov & Prodanc^ 
( l2005l) found a pseudo-potential for circular motion of a weakly 
charged particle in the Ker r-Newman space- ti me. A nother ex¬ 
tension was suggested by IStuchlfk & KovaH (l2008h who de¬ 
rived a generalization of the Paczyhski-Wiita prescription for the 
Schwarzschild-de Sitter black hole. 


In order to properly involve rotational d ragging, velocity- 
depen dent potentials have also been considered. ISemerak & Karas! 
( Il999l) tested one such idea against the exact solution on long-term 
behaviour of the di fference between t he respective free-particle 
dynamics. Recently iGhosh et alj ( l2014h suggested a new pseudo¬ 
potential which reasonably reproduces the Kerr space-time fea¬ 
tures for moderate centre’s angular momentum and moderate en¬ 
ergy of the orbiter (see also the overview given in Introduction of 
that paper, including previous results of its authors). But even in 
the Schwarzschild case the difference between Newtonian and rel¬ 
ativistic dynamics suggests the usage of velocity-dependent expres¬ 
sions; in a thorough stud y of the pseudo-Newtonian descriptions of 
the Schwarzschild field. iTeie da & RosswogI (l 2013l) brou ght such a 
more advanced possibility (see lTeieda & Ross wodl2014l for its fur¬ 
ther development). 
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values of i/M chosen (as going from bottom to top curves) 
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Figure 1. Comparison of effective potentials resulting from the pseudo-Newtonian gravitational potentials 0,0 and O (enlarged by one) with the exact 
Schwarzschild effective potential — 2M/r)(l -|- (?■ . Several profiles with different I are plotted, with the values of I adjusted (differently for 

different potentials) so that the curves be similar; particular curves involving the marginally stable and marginally bound circular orbits are shown for all 
the potentials and are easily recognizable. (For the NW potential, the £ = 0 case does not differ much from iJM = 1, so it is not included.) All the three 
potentials look similai' and not far from the actual Schwarzschild one. Our logarithmic potential is clearly very close to the PW one, having its maxima at 
slightly larger radii; the NW-potential profiles, on the contrary, are shifted to smaller radii with respect to the PW ones. The main difference is that the slopes 
of Schwarzschild curves are less steep. Also notice the differences in the values of I corresponding to roughly same heights of the potentials; i) those required 
for the NW potential are considerably lower; ii) at the high-i end, those required by the Schwarzschild potential grow faster. 


2 BLACK HOLE WITH DISC OR RING: A 
PSEUDO-NEWTONIAN DESCRIPTION 

Exact superpositions of a vacuum static axisymmetric (originally 
Schwarzschild) black hole with a concentric thin disc or ring are 
described by formulas which were given in the previous papers of 
this series (see mainly section 1.1 of the last paper III for a compact 
summary), so rather than repeating them again, let us only specify 
that we will again choose the inverted 1st member of the counter¬ 
rotating Morgan-Morgan thin-disc family (iMMl disc) and the 
Bach-Weyl linear ring (BW ring) as the external sources, approx¬ 
imating a thin accretion disc or ring, respectively. Let us also re¬ 
mind that (t, p, z, (j)) stand for the Weyl coordinates and {t, r, 0, cp) 
for the Schwarzschild-type coordinates, with t and 4> being Killing 
time and azimuth and p, z oi r, 8 covering the meridional two- 
space. Geometrized units are used in which c = G = 1, cosmolog¬ 
ical constant is (necessarily) set to zero and index-posed commas 
mean partial differentiation. 

Newtonian analogue of the relativistic black-hole-disc/ring 
picture studied in previous papers is given by function v which 
determines the gtt metric component, in Weyl coordinates satisfies 
the Laplace equation and represents a direct counter-part of New¬ 
ton’s gravitational potential. We will thus use the metric functions 
and zzbw of the disc and of the ring directly as the disc or 
ring Newtonian potentials, respectively. The Schwarzschild-centre 


potential izschw, on the other hand, will be just emulated by a cer¬ 
tain effective pseudo-potential. We will test three simple cases. 


Vpw = 

M 


(1) 

r - 

-2M 


= 

M 

0 - 

3M , 

(2) 

r 

+ 9 ) ! 

V / 

Uln = 




(3) 


The first was proposed by Paczvnski & Wiital ( ll98(]h . the second by 
iNowak & Wagon^ (Il991 ), and the logarithmic one represents an¬ 
other possibility we are submitting for comparison. The Paczynski- 
Wiita potential is a default benchmark, very simple yet behaving 
surprisingly well in many situations. The logarithmic potential is 
simply included because we newly suggest it here. And the Nowak- 
Wagoner potential is chosen for it has yet another form which will 
be seen to result in a rather different character of the accessible 
phase-space region; at the same time, it has turne d out to be the 
best o f “simple” possibilities in some studies (e.g. ICrispino et al.l 

i2onh . 


Other major simple pse udo-Newtonian subst itutes for 
Schwarzschi ld wer e provid e d by Artemovaet^ 1 1996h and quite 
recently by IWegd ( |2012|) . lArtemova et al.l ( Il996h used several 
pseudo-potentials in studying disc accretion onto black holes; in the 
non-rotating case, they considered expressions (we number them 
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according to the original paper) 


VabN3 ~ 1 + Y ^- - ’ 

Vabn4 = 5 In (l - ■ (5) 

The second one (just equal to the Schwarzschild potential t^schw) 
is similar in form to our logarithmic expression but we will see 
that the latter is actually more similar to the PW potential (see Figs. 
[2and[^. A comparison of t he two ABN potentials with the PW and 
NW ones was performed bv ICrisnino et al.l (120111) on the motion of 
a particle emitting scalar radiat ion. More re cently, several serious 
options have been presented bv lWeed )l2012h (original marking by 
A, B and C is kept again), 


Vwa = 

M 

r 

(-“)• 

(6) 

VwB = 

M 

r 

/ 3r 4M\ 

\3r — 5M 3r / ’ 

(7) 

Vwc = 

M 

1 + ^(3 -V6) + 2^(5 - 2v^) 

(8) 

r 

1 

1 


and shown to yield better results for the apsidal precession of low- 
energy (about parabolic) orbits than the Paczyhski-Wiita potential. 
Recently we have included, with a surprisingly good result, Vwa 
in a c omparison of li ght-ray approximations in the Schwarzschild 
field (is emerald l2015h . However, this potential turns out to be un¬ 
suited for our present purposes as shown in the next section (equa¬ 
tion (lH} and below). All the other four potentials, Vabns, Vabn4, 
VwB and Vwc, are included in Figs. and [3 in order to at least 
illustrate their basic nature against those we are going to study in 
more detail in the present paper. 

2.1 Issue of comparison in coordinates 

When preparing to superpose the centre-describing potentials with 
I'iMMi or vbw, one encounters the main query, however: How 
exactly to perform the Newtonian superposition in order to get a 
plausible counter-part of the relativistic system? Which coordinates 
covering the curved relativistic space-time are adequate counter¬ 
parts of Euclidean coordinates of the Newtonian description? The 
Newtonian pseudo-potential for the black hole is usually given in 
Euclidean spherical coordinates and simulates the hole represented 
in Schwarzschild coordinates, while the disc/ring potentials are nat¬ 
urally taken over from relativity in cylindrical coordinates. In the 
relativistic description, the linear superposition holds in Weyl coor¬ 
dinates p, z which are of cylindrical type and where the black-hole 
horizon appears as a finite line singularity at p = 0, | 2 | < M. Af¬ 
ter transformation to Schwarzschild-like coordinates of spheroidal 
type, 

p = \Jr{r — 2M) sin 9, z = {r — M) cos 9, (9) 

the black-hole horizon becomes spherical, while the disc/ring keeps 
its shape but has a slightly bigger coordinate radius. 

The spheroidal character of the black hole is clearly not well 
represented in the Weyl coordinates. However, since the relativis¬ 
tic superposition is performed in them, one should probably repro¬ 
duce it in the Newtonian approach in the following way: i) take the 
pseudo-Newtonian potential (in spherical coordinates) and trans¬ 
form it into the Weyl coordinates; ii) add the disc or ring potential 
expressed in the Weyl coordinates; iii) transform the result to the 



Figure 2. Values of the angular momentum i needed to raise the centrifu¬ 
gal barrier to a given energetic level £ (thus to establish an unstable cir¬ 
cular orbit with that energy), plotted for the potentials we compare {£ is 
enlarged by 1 for the potentials in order to match the relativistic case). The 
Nowak-Wagoner potential yields the worst result and our logarithmic po¬ 
tential yields the best one, yet none of them reproduces the Schwarzschild- 
field behaviour properly. The curves provided by potentials (D and m of 
Artemova et al. (1996) are also shown in dashed grey and the potentials m 
and i) of Wegg (2012) are drawn in dotted brown. 


Schwarzschild-type coordinates. Since the Newtonian fields super¬ 
pose linearly in any coordinates, one can summarize this without 
the intermediate step: take the black-hole pseudo-potential and add 
to it the disc or ring potential transformed from cylindrical to spher¬ 
ical/spheroidal coordinates in a Weyl-like manner, i.e. by substitut¬ 
ing (|9}. 

Alternatively, rather than to take over the transformation be¬ 
tween the Weyl and Schwarzschild coordinates to the Newtonian 
description, one could assume that the relativistic disc/ring poten¬ 
tial in Weyl coordinates corresponds to the Newtonian one in com¬ 
mon cylindrical coordinates, connected with the spherical ones (in 
which the pseudo-potential for the centre is given) by the Euclidean 
relation p = rsinS, ^ = rcos9. However, since the pseudo¬ 
potentials should imitate the black hole, which means mainly to im¬ 
itate the occurrence of the horizon, it is reasonable to demand that 
the spheroidal-cylindrical transformation have in both cases simi¬ 
lar effect on the central source: if it shrinks the relativistic source 
into a rod, it should not leave the Newtonian source intact (as the 
Euclidean-type relation). We have anyway tested this second pos¬ 
sibility too and learned that if the external source is not very close 
to the centre (below lOM, say), the results are almost identical. 

However, carefully as one may try to consider the correspon¬ 
dence between the relativistic and pseudo-Newtonian systems, they 
inevitably remain different, the more so that not only the space(- 
time) backgrounds differ, but also the dynamics (equations of mo¬ 
tion), so one should at least expect a quantitative discrepancy, un¬ 
less employing some more sophisticated velocity-dependent poten¬ 
tial. 
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3 MOTION IN MODIFIED NEWTONIAN POTENTIALS 


The motion of test particles in the velocity-independent axially 
symmetric Newtonian potential V (r, 9) is described, in spherical 
coordinates (r, 9, (j)) and with obvious notation, by equations 

f = —Vr + r [9"^ + sin^ 9), (10) 

9 = sin 9 cos9 — 2rr9 — Vg , (11) 

<j> =—2rij> [f + r9 cot 9). (12) 


If the field is even spherically symmetric, V = V{r), the term 
in the 2nd equation vanishes and the motion gets confined to a 
plane. The orbital plane is usually identified with 9 = -k 12, so one 
is left with equations 

r = —V,r + r(j>'^, r(j> = —2^r. (13) 


These have energy and angular-momentum integrals 

i? = ^ (f^ + + mV, L = mr^ ij>, (14) 

which invert for velocities as 


I .2 2mr\E-mV)-L^ , ,,,, 

^2’ ^ = - ^2 - = 2(i:-Veft), (15) 


wher^il 


Vefr:=V^+:^, £■■=-, (16) 

2r^ m m 


Circular orbits exist where 


Veff,r “0 <=> 7 


so their linear speed amounts to 

rb = , 


(17) 


(18) 


their energy is given by the corresponding potential value 

£{f = r^V,r) = Vefi(^" = r-®V;.) = V + ^rV,r (19) 
and their stability is determined by the sign of 

OT/ 

V,S,rr{f = rW,r) = V,rr+^^. ( 20 ) 

r 

The character of radial motion and its response to perturba¬ 
tions are thus governed by shape of the potential well (given by V 
and 1) and by the particle’s specific energy £. Most importantly, 
the shape of Vefi and the value of £ determine the properties of the 
region accessible to the particle within the (r, f) diagram. A well 
known crucial point is whether this region is closed or open towards 
the centre, which, for a given energy, depends on the height of the 
centrifugal barrier id jr^. In the marginally closed state, the acces¬ 
sible domain is bounded by a separatrix which corresponds to a ho¬ 
moclinic orbit, winding - in infinite past and infinite future - from 
and on the unstable circular periodic orbit residing at the potential 
saddle-point vertex. Homoclinic orbits, a salient feature of black- 
hole fields, represent an infinite-whirl limit of the zoom-whirl type 
of motion (a strong-field bound motion with extreme pericentre ad¬ 
vance), and are familiar to mark the frontiers of chaotic regime - 
their perturbation leads to the occurrence of a “homoclinic tangle”, 
through which the original circular orbit breaks up into a fractal set 
of periodic orbits. 


^ As noted in figures and their captions, we actually shift the specific en¬ 
ergy £ by one so that a particle at rest at infinity has £ = 1 in accord with 
the relativistic case. 



_ i.Nir. _^_._I_ Zj _ Zj _^_._I_._._._^_ 

0 10 20 30 


Figure 3. Top: One specific effective-potential profile plotted for all the 
gravitational potentials considered, with the angular momentum I = 
3.75M (this value is chosen in most of the figures presented in next sec¬ 
tions). Like in Fig. |2] the Artemova-et-al. potentials Vabn 3 and Vabn 3 
are also shown in dashed grey and the Wegg potentials VwB and Vyvc 
are drawn in dotted brown. Only the PW and the In potentials (red and 
green) seem to approximate the exact Schwarzschild shape in some way. 
Clearly the PW potential is more open towards the centre, while the In po¬ 
tential is more closed than the actual Schwarzschild case. Bottom: Similar 
plot, but with the angular-momentum values chosen so that all the effective 
potentials have the same maximum £ + 1 = 0.987746 at the unstable cir¬ 
cular orbit (in the Schwarzschild case, one takes just £). Concretely, this 
means £ = 3.9M for Schwarzschild, £ = 3.9494M for Paczynski-Wiita, 
£ = 2.7475M for Nowak-Wagoner, i = 3.7805M for the logarithmic 
potential, £3 = 2.5739M and £4 = 3.1028M for the Artemova-et-al. po¬ 
tentials Vabn 3 and Vabns, and £b = 3.9651M and £c = 3.9735M for 
the Wegg potentials Vwb and Vwc ■ The pseudo-potentials yield somewhat 
different radii of the unstable circular' orbit (only the PW and In potentials 
have it very close to the con'ect value) and their valley existing above this 
orbit is deeper than the Schwarzschild one; the difference is especially large 
for the ABN potentials. 


The homoclinic orbit is infinite, but the length of its trail in 
reasonable coordinates (r, r in our case), i.e. of the accessible- 
region bounding separatrix, indicates the size of a phase-space re¬ 
gion which turns chaotic under perturbation. This does not provide 
any plausible (“covariant”) measure of what fraction of the phase 
space will be affected, but still can be used to compare different po¬ 
tentials. A similar suggestion (only given by d rather than by f) is 
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contained in the length of the potential valley Veft(^circ) below the 
energy level of the unstable circular orbit or in this valley’s area. 

Let us now briefly check the basic properties of effective po¬ 
tentials given by the gravitational potentials in particular 

whether and how they reproduce circular periodic orbits, decisive 
for the response of the dynamical system to perturbation. However, 
consider first the Wegg’s expression lUl in order to realize why it is 
not suitable this time. For the corresponding effective potential. 


Veft = — 


M 

r 




( 21 ) 


the condition for circular orbits P = r^V,r yields Mr — P — 
6M^ for the radius, so P > 6M^ must hold in order that such 
radii really exist. But for the Wegg potential one has Veft.rr = 
{P — 6M'^)jr*, so all the P > 6M^ circular orbits sit at the 
potential minimum, hence they are all stable and not interesting for 
us. Therefore, rather than mimicking the occurrence of unstable 
circular orbits, so characteristic to the black-hole fields, the Wegg’s 
A-potential behaves like Newtonian —M/r, just with the critical 
value of P shifted from zero to 6M^. (This is no wonder, since 
Wegg suggested the potential specifically for near-parabolic orbits 
at larger radii.) 

The shapes of the effective potentials resulting from the 
Paczynski-Wiita, Nowak-Wagoner and our logarithmic potentials 
is compared in Fig.[T] All the three potentials host both stable and 
unstable circular orbits and are clearly quite similar. They all yield 
the correct radius r = QM for the marginally stable circular or¬ 
bit (ISCO). The Paczyhski-Wiita potential also does so for the 
marginally bound orbit (IBCO, r = AM), reproducing besides the 
angular-momentum Schwarzschild value i = AM there. On the 
other hand, the logarithmic potential gives the correct value of an¬ 
gular momentum at the ISCO {t = 2v/3M). The latter is a conse¬ 
quence of a more general tuning: circular orbits of the logarithmic 
potential satisfy 




Mr'^ 
r — SM 


( 22 ) 


which is exactly the same expression as would be obtained in the 
exact Schwarzschild field. This means, in particular, that a Keple- 
rian disc in the In-potential would have exactly the same distribu¬ 
tion of angular momentum as in the Schwarzschild case. 

Figure emphasizes what may not be evident from Fig. [T] 
that although the shapes of the potentials seem similar to the actual 
Schwarzschild one, they may differ significantly or just fail in some 
important aspects like the relations between the energy and angular 
momentum for the unstable circular periodic orbits. Specifically, a 
particle with £, I located below the respective curve in Fig. will 
orbit in an allowed region open towards the center and will thus be 
prone to black hole in-fall; on the other hand, particles from above 
the curve will orbit in two distinct regions, the exterior one being 
closed-off from the center by the centrifugal barrier. However, if 
one picks £{+1) < 1 (hoping for bound motion later harbour¬ 
ing chaos) and £ too far above the curve, there might be no bound 
particles orbiting the black hole because of a too high centrifugal 
barrier. Hence, in the f (-1-1) < 1 range the PW, Wegg B and C, 
and log potentials are expected to exhibit satisfactory behaviour in 
terms of the overall nature of the allowed region, whereas the NW 
and Artemova potentials will not show a good correspondence. 

One can judge from this that although the character of chaos 
induced by perturbation of the pseudo-Schwarzschild fields is 
likely to be similar to what is a common experience from weakly 
non-integrable systems, its dependence on the relevant parameters 


will be quantitatively different, in particular the parameter values 
critical for an occurrence of various features (resonances, separatri- 
ces, chaotic layers) will be different. Also, as the potential valleys 
provided by the pseudo-potentials are generically deeper than the 
actual Schwarzschild ones (see Fig. 13, it might be loosely antici¬ 
pated that the corresponding Newtonian motion will rather be more 
chaotic than geodesic motion in the exact relativistic field. How¬ 
ever, one must remember that we are yet talking about the central 
black hole only, and, also, that the relativistic dynamics is different 
from Newtonian (already special-relativity effects make some dif¬ 
ference), so the centre’s effective-potential shape is just one part of 
the story. 


3.1 Superposition of the black hole with a disc or ring 


The second part is the gravitational potential of the additional 
source which in our case will be represented by a thin annular disc 
or a ring. If a static and axially symmetric source is placed around 
the centre, the field is no longer spherically symmetric, hence a 
generic motion is no longer plane-like and one must return to equa¬ 
tions Gll-Glll. Their energy and angular-momentum integrals now 
have the form 

£ — — {'P A- r^(P A- sin^ 0) A-V, £ = r'^ (j) sin^ 6 , (23) 
and invert for velocities as 

0= =2(£--VefF), (24) 

r-^ sin 0 

where 


Veff := y -b ■ 


2r2 sin^ 0 ' 

To obtain an effective potential for the motion in the complete 
field of the (pseudo) black hole surrounded by some external source 
(which generates potential iZext), one simply takes the above Veu 
with 

V = Fpseudo(r) A- Oext ^a/ r{r — 2M) sin 0, (r — M) cos . 

We illustrate the possible outcome by adding the inverted 1st 
Morgan-Morgan counter-rotating disc which was already involved 
in previous papers of this series and whose gravitational potential 
reads, in the Weyl-type cylindrical coordinates 0, 


^disc — 


M 


7r(p2 + 22)3/2 


n 2 . r, 2 ,2 P ~2z 

+2^ 

2 


-i(3E-36"-bp"-b2") 


arccot S 


(26) 


(see e.g. lZacek & Semet^l2002h . where 

E sj (p2 — 6^-1- 2^)2 A- Ab'^z'^ 


S ■- 


E - (p2 - + ^2) 

2 (p2 -b 2;2) 


and AA and b denote mass and Weyl inner radius of the disc. Figure 
|4] shows the results obtained with different pseudo-potentials and 
compares them with the one following from an exact relativistic 
treatment which describes the problem (geodesic motion) by equa¬ 
tions (see section 4 of paper I) 


g 2 (A Ag^hw) -I- r{r — 2M)(M®)^j = £^ — 


(Vefi)^ ;= 




£ 2 g 2 i/,ji 3 „ \ 

sin^ 0 ) 


g^^disc 


(Veft)", (27) 
(28) 
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Figure 4. Meridional (c/i = const) sections of the effective potentials for an originally Schwarzschild black hole surrounded by the 1st member of the 
Morgan-Morgan counter-rotating thin-disc family. Exact relativistic superposition is shown (1st row) together with those involving Paczyhski-Wiita (2nd 
row), logarithmic (3th row) and Nowak-Wagoner (4rd row) imitations of the black hole. All the cases are determined by the value of specific energy £{-|-l) = 
0.987746 at the unstable circular orbit as in Fig.[^(so the values of I are also exactly as there). Within all of the four rows, the disc relative mass MjM is 
chosen, from left to right, 0.0, 1.0 and 5.0. In all the plots, the contours shown are Vefj(+1) = 0, 0.1, 0.2, 0.3, ... 0.700, 0.705, 0.710, 0.715, ..., 0.990, 
0.995, 1. (Vefj -b 1 is taken in the Newtonian cases, whereas just Vgfj in the Schwarzschild one.) Axes are scaled in the units of M. 


where A has to be computed by line integration of the gradient of 
total potential iz, with Aschw denoting its pure-Schwarzschild form, 
is four-velocity of the test particle, and £ := —ut and £ := 
are constants of the geodesic motion following from the Killing 
symmetries (they represent specific energy and specific azimuthal 
angular momentum of a test particle with respect to infinity). The 
figure confirms that the pseudo-potentials we consider here pro¬ 
vide similar - but not the same - effective potentials as the exact 
Schwarzschild-field description, with the Paczynski-Wiita and our 
logarithmic formulas apparently being quite close to each other. 

Superposition with the Bach-Weyl ring is acquired in the 
same manner, just with tZgxt represented by 


where K{k) ~ ^ ■ 2 complete Legendre el¬ 

liptic integral of the 1st kind and A4 and b are mass and Weyl radius 
of the ring. 


4 NUMERICS 

Trying to check our previous results also by using a different nu¬ 
merical method(s), we turned to symplectic integrators, suitable for 
conservative systems. However, the two outer sources we consider 
differ in what to do when the particle hits them: the ring is a cur¬ 
vature singularity, so it is appropriate to halt the trajectory if it gets 
to its closest vicinity, whereas the thin disc is only singular at its 
inner edge while cross transitions elsewhere are approximated as 
non-collisional (pure gravitational effect). Hence, the disc case has 
to be treated more carefully, regarding that there is a normal-field 
jump across the equatorial plane (hence jump of the z component 
of acceleration) above the disc inner radius. 

More specifically, the geodesics in the fields given by superpo¬ 
sitions with the Bach-Weyl ring are integrated using the 6th-order 
explicit symplectic partitione d Rung e-Kutta method with coeffi¬ 
cients adopted from lYoshidal ( 1 19901) (Solution A) and with step 
= (2 ^ 5) ■ 10“^M depending on the strength of the ring. 

In the case of thin discs (1st inverted Morgan-Morgan disc 
in our case here), regular integrators bring linear to polynomial 
growth of error in constants of the motion due to the jump in 
vertical field. In previous papers of this series, we got over this 


izbw = — 


2M 


K 


TV \^/ip+b) 


2 
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by the Hut’a method with adaptive step and using higher float 
representation. For the present paper, we developed a different 
variable-step integrator largely inspir ed by the IGEM code of 
Sevrich & Lukes-GerakopoulosI |20l3) and having the desirable 
properties of reversibility and symmetry (see lStofferl 19951 for other 
varible-step symmetric-reversible integration methods). It is based 
on Gauss collocation method with 3 collocation points (s = 3) and 
step size determined by the collocation points. Unlike in IGEM, the 
step size is not determined by the Jacobian of the integrated vector 
field but by spatial coordinates and the integrated vector field it¬ 
self (i.e. by phase-space variables and their time-derivatives at the 
collocation points). 

We start by choosing the step 


ho 


e 

11/(1) + /(s)ll 


(30) 


where f ~ f{xi) is the integrated vector field at points Xi 
of the Gaussian collocation and the norm is defined by ||/|| := 
X] 1/^ |, where are components of the vector /. (Any norm actu¬ 
ally works. We use absolute value which is computationally less de¬ 
manding than fractional powers, for example). The integrator will 
be reversible and symmetric if h{xi) = h(p^,... ... ,q^) 

is a function symmetric with respect to the reversal of order, 1 -f-)- s, 
and to the change of the sim of momenta, p^ —> —p^■ Now, the step 
ho is adapted according tqf| 


where 


ho _ e 

«(C) «-(C) 11/(1)+/(s) 


n(C) := 1 -f 


^1 

52+e' 


c 



(31) 


(32) 


so it remains about ho for tii, while for ^2 + < ^i it 

is contracted; the factor which multiplies ho is however never less 
than (52/5i. The coefficients 5i,S2 are set so that the particle travels 
in a controlled manner as close to the equatorial plane as possible. 

Then, from some minimal C,, the particle is reflected with re¬ 
spect to the equatorial plane: when its |(^| falls below some chosen 
Cmin, the program first estimates whether it will cross the equatorial 
plane in the next k steps by computing 

thus basically using the Euler explicit method with a step of roughly 
K,ho \ if C is found to change sign, the original position is reflected 
by ^ —C,. The advantage of this approach is that the particle 
encounters a “stepping wall” near = 0, the iterative Gaussian 
collocation does not suffer from the nearby discontinuity and the 
— C reflection exactly conserves energy. The only point vio¬ 
lating the integrator’s symmetry is the step estimate of the cross¬ 
ing, but any symmetric reversible stepping would be implicit and 
difficult to iterate over the discontinuity, with only small bene¬ 
fit to accuracy. We checked that when the parameters are tuned 
properly, the error typically oscillates without any drift, as typical 
for symplectic/reversible-symmetric integrators. In some cases the 
self-adjustment of the step has proven insufficient and a slow lin¬ 
ear growth in relative energy error was observed (usually for parti¬ 
cles infalling onto a black hole), but this error only rarely exceeded 


^ We perform the integration in Euclidean r sin 6, r cos 6 (not in the Weyl- 
type coordinates), so we better introduce C = r cos 9 z) to avoid con¬ 
fusion. 


10 By numerical experiments, we have found the following 
parameter ranges to be optimal: 

Cmin = (1 A 5) • 10-^M, 

51 = (10"®^ 10”®) 

52 = (10"®^ 10“®) 
e= (5 3-8) • 

«;= 13-3. 

Let us add that the Gaussian collocation was found by fixed- 
point iteration and convergence was confirmed by checking the dif¬ 
ference between the current set of collocation points Xi and the 
previous one as represented by 

s 2JV 

^ = EEhi-4p)|’ (34) 

i=l j = l 

where W = 2 is the number of degrees of freedom. The iteration 
was stopped whenever A < 10“^®. Such a tolerance corresponds 
to an average error of the order of 10”'^'^ per collocation compo¬ 
nent, which is about what can practically be achieved, because spa¬ 
tial position (configuration part of x) was often larger than 10, the 
“distance” A includes subtraction of close numbers and we used 
double precision which stores about 15 digits. 

The Poincare surfaces of section where created of 3600 
equatorial-plane intersections, recording transitions in both direc¬ 
tions. Whenever the singularities of the central potential or the ring 
were closely approached, the integration was stopped and restarted 
again with a nearby initial condition until a sufficient number of 
points was collected. However, the whole set of intersections gener¬ 
ated by a given trajectory was discarded if a relative error in energy 
turned out to be too large (namely > 10“®). Overall, the initial 
conditions were chosen by a pseudo-random algorithm similar to 
the one described in paper I. 


5 COMPARING EXACT RELATIVISTIC AND 
PSEUDO-NEWTONIAN PICTURE 

Let us stress once more that the (pseudo-)Newtonian and relativis¬ 
tic dynamical systems in question are fundamentally different, be¬ 
cause they live in a different configuration space and their evolution 
is described by a different dynamics as well. It is even impossible 
to decide which situations are “similar”, because most of the rel¬ 
evant variables actually have different meaning within Newtonian 
and relativistic case; for example, if one places the external source 
at some “given” radius, it has different meanings in the Euclidean 
spherical/cylindrical coordinates r sin 9, r cos 9, in the Weyl-type 
coordinates p = y/r{r — 2M) sm9, z = {r — M) cos 9 (which 
we use here) and - in relativity - in terms of proper radial distance 
or in circumferential radius. Therefore, one can only wish for a rea¬ 
sonable correspondence of qualitative phase-space features and of 
their evolution with analogous parameters. Yet it will still be inter¬ 
esting to see whether and which of the potentials reproduce at least 
some of the quantitative aspects, like the pattern of resonances and 
the sequence of their appearance. 

Needless to say, one has only a very restricted space here for 
such a comparison. It is symptomatic for non-integrable systems 
that their dependence on parameters is “chaotic” (non-smooth) - 
they may change only slowly within one parameter range, whereas 
very abruptly within the other (which may be quite narrow). Being 
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Figure 5. Poincare diagrams in axes (r, v'^) showing passages of geodesic orbits with conserved energy S + 1 = 0.955 and angular momentum i = 3.75M 
through the equatorial plane of a centre described by the Paczyhski-Wiita potential (with mass M) and surrounded by an iMMl disc with inner radius 
r^isc = 20M. Dependence on mass of the disc A4 is shown, as given in the plots. Accessible sector is indicated in purple and r axis is in units of M. 
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5 10 15 20 25 10 20 30 40 50 


Figure 6. Same series of plots as in Fig. [5] but with the central black hole simulated by the logarithmic potential Comparison of these two figures with 
fig. 4 of paper I indicates that the phase-space portrait of all the three dynamical systems is similar, though various quantitative differences can be noticed (see 
mainly behaviour of the accessible region) as more discussed in the main text. 
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Figure 7. Poincare (r, v'^) diagrams showing passages of geodesics with angular momentum £ = 3.75M through the equatorial plane again, for the centre 
described by the Paczyhski-Wiita potential (with mass M) and suri'ounded by an iMMl disc with mass A4 = 0.5M and inner radius rdisc = 20M. Here 
dependence on energy of the orbiters is in focus, as indicated in the plots (we enlarge it by unity to match with the relativistic value). 
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Figure 8. Same series of plots as in Fig. [7] but with the central black hole simulated by the logarithmic potential m- Compaiison of these two figures with fig. 
6 of paper I again indicates that both Newtonian dynamical systems well approximate the relativistic one; quantitative differences are further discussed in the 
main text. (Mainly evident is the different delimitation of accessible phase-space sector again, following from differences in effective-potential profiles.) 
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Figure 9. Poincare diagrams in axes (r,v’~) showing passages of geodesic orbits with conserved angular momentum £ = 3.75M through the equatorial 
plane of a centre described by the Nowak-Wagoner potential (with mass M) and surrounded by an iMMl disc with inner radius ruisc = 20M. The first two 
rows show dependence on mass of the disc A4, while all the orbiting particles have energy £ + 1 = 0.955. The last row shows just 3 examples of how the 
plots change with orbital energy £, while the disc mass is set at M = 0.5M. The plots are rather different from those involving the Paczyhski-Wiita or the 
logarithmic potential, because the Nowak-Wagoner potential is so weak that it does not form “its own” valley and the accessible region is maintained by the 
disc, at least in the Ad-series plots. On the other hand, exactly due to this different character it is useful to include at least this one series employing the NW 
potential. 


only able to select several sections through the very rich parame¬ 
ter space of the systems, one can either take those with the same 
values of the corresponding parameters, those showing similar fea¬ 
tures, or simply those where something interesting is happening. 
Without adhering to any strict rule, we have generally set the fixed 
parameters at formally identical values as in the relativistic case 
and varied the free parameter in roughly the same range. 

The comparison in general reveals that the overall tendency is 
the same in both the relativistic and pseudo-Newtonian systems: 
when the perturbation strength (disc mass in our case) or parti¬ 
cle energy increases, the system first gets more and more chaotic, 
whereas for very large parameter values the “primary” regular re¬ 
gion slowly grows again. However, since such a behaviour is quite 
typical for weakly non-integrable systems, we will mainly try to 
note the differences. 

We start by evolution of the phase portrait with mass of the 
external source. Figures and show how Poincare diagram of 
equatorial transitions changes with relative mass of the inverted 
1st Morgan-Morgan disc as the external source. Placing the in¬ 
ner edge of the disc on raise = 20M and setting £ = 0.045 (= 
0.955 — 1.000), I — 3.75M as in paper I, the figures present di¬ 
agrams obtained for 8 different values of MjM between O.IM 


and 1.7M. The Schwarzschild centre is imitated by the Paezynski- 
Wiita potential in Fig. Awhile by the logarithmic potential in Fig.|^ 
We have not included the Nowak-Wagoner potential in the detailed 
study, because it has turned out to yield rather different results, not 
well compatible with the exact relativistic picture (the NW poten¬ 
tial is “too weak” and for a large portion of the studied parameter 
ranges its phase space bears no bound particles); however, Fig.j^is 
provided for cursory illustration. 

The Figs.[3and|^are to be compared with fig. 4 of paper I. 

• The main difference concerns the accessible domain which, in 
comparison with the exact relativistic case, is more open towards 
the centre for the Paezynski-Wiita potential, whereas more closed 
for the logarithmic potential (see Fig. [3: for real Schwarzschild, 
the domain is closed first, enlarges with increasing disc mass and 
finally opens towards the centre when the disc mass reaches about 
half of the black-hole mass (this applies specifically to the iMMl 
disc with raise = 20M, of course). In contrast, for the Paezynski- 
Wiita potential the accessible sector is always open towards the 
centre, whereas for the logarithmic potential it is closed and only 
opens after the disc outweighs the centre. However, this does not 
seem to be that crucial for evolution of the phase-space features, 
the opening only enables the centre to “suck out” the outer chaotic 
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Figure 10. Poincare (r, v'^) diagrams showing passages of geodesics with conserved energy S- -\- 1 = 0.977 and angular momentum i = 3.7bM through 
the equatorial plane of a centre described by the Paczyhski-Wiita potential (with mass M) and surrounded by a Bach-Weyl ring with radius rring = 20M. 
Dependence on mass of the ring is shown, with values given in the plots. Accessible sector is indicated in purple and r axis is in units of M. 
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Figure 11. Same series of plots as in Fig. IlQI but with the central black hole simulated by the logarithmic potential (^. Comparison of these two figures with 
fig. 10 of paper I indicates that the phase-space portrait of all the three dynamical systems is similai; though many quantitative details are different (it would 
be pointless to discuss them extensively due to the richness of the structure); note again the different delimitation of the accessible region. 
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Figure 12. Poincare (r, n’’) diagrams showing passages of geodesics with angular momentum (. = 3.75M through the equatorial plane again, for the centre 
described by the Paczyhski-Wiita potential (with mass M) and surrounded by a Bach-Weyl ring with mass M = 0.5M and radius r^ing = 20M. Here 
dependence on energy of the orbiters £ is in focus, as indicated by its values given in the plots (we enlarge it by unity again). We are not showing plots obtained 
for £ + 1 = 0.910 and less which only contain a tiny accessible region around the ring (the other region between the centre and the ring is not existing yet). 


sea. (This makes the open diagrams asymmetric with respect to 
w" = 0.) 

• The similarity of all three systems is really striking, with most 
phase-space structures appearing and in the same succession. In the 
Paczynski-Wiita case, similar features appear for somewhat lower 
disc-mass values (about 0.1M^0.2M “earlier”) than in the rela¬ 
tivistic case, while for the logarithmic potential they appear still 
about O.OSM^O. IM earlier than for the PW potential. This may be 
interpreted as slightly stronger inclination of the pseudo-Newtonian 
system towards chaos, which is in accord with our preliminary 
guess stemming from deeper potential valleys provided by them 
(see Fig.l^l. 

• More details about the structures: with increasing perturbing 
mass, the relativistic geodesic system (see paper I) first develops 
a 3-fold island within the primary regular region (2:3 resonance, 
“fish”-shaped orbit in Fig. 1171: then (temporarily) a 4-fold one ap¬ 


pears within the chaotic periphery of the accessible region: this is 
a particularly shaped “symmetrized set” of 1:2 resonances (anal¬ 
ogous feature appears “earlier” in the PW case)0 Later the cen¬ 
tral regular region gives birth to five islands (4:5 resonance, again 
identical in the relativistic and pseudo-Newtonian case), then even 
7-fold and 9-fold “baby-islands” (6:7 and 8:9 resonances) can be 
spotted, and finally the region breaks up into two parts symmetrical 
with respect to u” = 0 which disappear shortly after the disc mass 
reaches about the black-hole mass. Meanwhile, a central regular 
sector appears and grows gradually with the disc mass increased 


® Normally, an m:k resonance is associated with a fc-fold (fc-periodic) is¬ 
land. It is not clear whether the 4-fold island represents a tangent or a pitch- 
fork bifurcation of the 1:2 resonance (cf. also the following commentary on 
a 1:1-resonance bifurcation). 
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Figure 13. A counterpart of Fig. ll2l showing the same series of dependence-on-energy plots for the centre described by the logarithmic potential s All the 
parameters are kept from previous figure, i.e. £ = 3.75M, M = 0.5M, r^ing = 20M, and also the values of E are chosen equally, as indicated in the plots. 
In addition, we have kept exactly the same axis ranges, so the two figures can be compared easily. Their relativistic counterpart is fig. 12 of paper I. 


still more. 

The Paczyhski-Wiita centre with the iMMl disc also first give birth 
to the 3-fold island and then to the 5-fold, 7-fold and 9-fold ones, 
corresponding to the same resonances as in the relativistic case; the 
4-fold structure only appears in a light-disc stage (along the bor¬ 
der of the regular domain). The logarithmic potential yields very 
similar behaviour, with the 4-fold structure not occurring at all. 

• The breakup of the original central island is a very character¬ 
istic feature of the relativistic as well as of the pseudo-Newtonian 
systems; in all cases it occurs when the disc mass Ai reaches about 
that of the central hole (M). More specifically (Fig. □1: if one takes 
any point r, 9, r, 9 on the original central orbit (red) and applies the 
reflection 9 ^ n — 9 and/or velocity-reversal r —>■ —f, 9 —>■ —9, 
the same central orbit is obtained, just in a different phase. Namely, 
the central orbit is - up to a phase shift - symmetric with respect to 
reflection and reversal which are discrete symmetries of the Hamil¬ 
tonian. However, this symmetry of the whole phase space need not 


be respected by individual invariant structures. The multiplication 
of resonant islands is then a kind of “spontaneous symmetry break¬ 
ing”, because as the central orbit shifts to the strongly perturbing 
disc edge, it looses stability and bifurcates into two (green) orbits 
which are reflection-symmetric when taken together as a “sym¬ 
metrized set” (the reflection operation maps the points of the first 
trajectory on the second one and vice versa). These green trajec¬ 
tories later bifurcate even further in the radial direction, into 2 - 1-2 
“reversible-asymmetric” trajectories (blue and purple). The four 
small islands in the Poincare diagrams with M = 1.7M in Figs. 

and 1^ thus correspond to a symmetrized set of four distinct 1:1 
resonances. 

• Let us point to one specific difference finally: In the log- 
potential system, one observes a strong 5-fold structure correspond¬ 
ing to a 4:5 resonance inside the central regular region, existing 
from Af = 0.33M to A1 = 0.62M (we mean the one oriented 
so that one “vertex” island lies on the = 0 axis and towards the 
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r 


Figure 14. A scheme of the 1:1-resonance bifurcations occuning both in 
the relativistic and pseudo-Newtonian hole-disc field. The left part indi¬ 
cates the imprints of the trajectories on the Poincare section, while the right 
part illustrates the corresponding trajectory shapes (coloured respectively) 
within the (r sin 0, r cos 6) plane, with the dotted line always representing 
the equatorial plane 0 = tt/2. Note that the sections of the blue and purple 
trajectories would in fact depend on the direction of velocity on the curve. 
The original 1:1 island first breaks up “vertically” and then both new is¬ 
lands further decouple “horizontally”. This typically happens in stage when 
the phase space is the most chaotic, which in terms of the disc mass as the 
perturbing agent means ~ M. In the present paper, the three phases can 
be seen in Figs. |5]and|6] in plots with M = 0.9M, l.OM and 1.7M; the 
corresponding relativistic situations were plotted in fig. 4 of paper I (see the 
plots with M = l.OM and I.IM there; the 3rd phase was not shown). 


centre); in the PW-potential system, the similar structure is weaker 
and only persists from A4 = 0.54M to Af = 0.67M; in the exact 
system it does not appear at all. (However, it can appear rarely for 
a different type of disc and/or for a disc placed on different radius - 
see fig. 5 in paper I, plots with raise = 14M and 15M.) Notice also 
how in the pseudo-Newtonian cases that structure finally switches 
over to a complementary/reverse 5-fold pattern, with one “vertex” 
lying away from the centre (which is common in the exact system). 

Now we proceed to energy which is one of the most important 
parameters of any dynamical system. Figures 0 and [8] show how 
Poincare diagram of equatorial transitions changes with conserved 
energy of the freely orbiting test particles. Placing the “iMMl” disc 
of mass M = 0.5M from raise = 20M and setting £ = 3.75M 
as in paper I again, the figures present diagrams obtained for 8 dif¬ 
ferent values of £ between f + 1 = 0.95 and f + 1 = 0.98. The 
Paezynski-Wiita potential is used in Fig. 0 while the logarithmic 
potential in Fig. [8] The Nowak-Wagoner potential is only illus¬ 
trated briefly in Fig.[^ 

The Figs. [7] and [8] are to be compared with fig. 6 of paper I. 
The latter shows less stages than we present here, but the compari¬ 
son anyway confirms what has already been observed above in fig¬ 
ures illustrating dependence on the perturbing mass (the sequences 
in fact resemble the previous ones): the pseudo-Newtonian systems 
well simulate the exact relativistic one, they are just slightly richer 
of tiny structures and display major features somewhat “earlier” in 
terms of the relevant parameter (here energy). In this sense, they 
can again be called “more chaotic” than the exact system, with the 
logarithmic potential perhaps being slightly more prone to irreg¬ 
ularity than the Paezynski-Wiita one. In the left column of Fig. 
[S] notice the nice (center-vertexed) 5-fold pattern and its switch¬ 
over to the “complementary” pattern between f + 1 = 0.957 and 
£ + 1 = 0.958. 

The same kind of illustration - dependence on external mass 
and on orbital energy - is also provided for the black-hole-like cen¬ 


tre surrounded by a Bach-Weyl ring (with radius rring = 20M). 
Figures Go} and GD show how the equatorial (r, n’’) section 
through the phase-space evolves with relative mass of the ring, 
while energy and angular momentum integrals are chosen f -F 1 = 
0.977 and £ = 3.75M; the centre is described by the PW po¬ 
tential in Fig. [T^ while by the In potential in Fig.[TT] Figures [T2l 
andll3lshow dependence on energy of the orbiting particles, while 
£ = 3.75M and the ring mass is set at Al = 0.5M; again the PW 
potential is employed in the first figure, while the In potential is 
employed in the second one. Figures Gil and GD are counterparts 
of fig. 10 in paper I, while Figs.ll2landll3lare counterparts of fig. 
12 in paper I. 

The comparison with paper I again verifies quite close simi¬ 
larity of all the three dynamical systems (the exact relativistic one 
and those with the black-hole simulated by the PW or the In poten¬ 
tial). One might however notice many unlike details, but they are 
not worth careful discussion, because in the ring case the dynamics 
is apparently very rich of tiny structures, both regular and chaotic. 
The rich ornamentation follows from close encounters with the sin¬ 
gular source, so in future work - whether within exact description 
or using pseudo-potentials - it will be sensible to rather consider 
orbits not closely interacting with the ring, i.e. to choose the acces¬ 
sible region so that not to involve the ring radius. Such a configura¬ 
tion will also be more realistic since there are no literally singular 
sources in nature (cf. paper III where this point was checked in 
simple modelling of Galactic circumnuclear rings). 

Anyway, comparison of Figs. 1 121 and [T3] with fig. 12 of pa¬ 
per I indicates, similarly as the centre-disc plots above, that the 
pseudo-Newtonian imitations of black hole lead to slightly faster 
“evolution” with parameters than the exact relativistic case, which 
can perhaps be interpreted as more “unstable” response to the per¬ 
turbation. For instance, the breakup of the principal regular sector 
existing below the ring occurs at f = 0.945F-0.950 in the relativis¬ 
tic case, while at £" +1 = 0.940 A 0.945 in both pseudo-Newtonian 
cases. Again quite different is the moment of opening of the ac¬ 
cessible domain towards the centre: in the relativistic picture, this 
happens at f ~ 0.934, while with the PW potential it happens at 
f + l~0.917 (the PW potential is almost ever open) and with the 
In potential it happens only at f + 1 ~ 0.953 (the In potential is 
almost ever closed). 


5.1 On dependence on perturbing mass and on orbital 
energy 

The Newtonian dynamics allows for a straightforward and quan¬ 
titative explanation of the correspondence between the changes in 
sections caused by variation of the perturbing mass M and by vari¬ 
ation of the orbital energy £. First, as we fix the total energy £ and 
increase the disc mass A4, the potential well becomes deeper, so 
the particle necessarily gets more kinetic energy. Hence, although 
the parameters of the surfaces of section in Figs. and might 
look like we study “identical” ensembles of trajectories subjected 
to a stronger and stronger dynamical perturbation, effectively it is 
not so. 

To illustrate this point further, we compute the average speed 
v(£,M,£, raise) over the equatorial plane for the parameters cho- 
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Figure 15. Average speed {35) with which the orbits (having I = 3.75M) 
intersect the equatorial plane of the system of a black hole surrounded by 
the iMMl disc with inner radius “rdisc = 20M. The dependence of v on 
the relative disc mass AijM is plotted for S + 1 = 0.955 (these curves 
are given in x crosses; top axis applies); the dependence of v on the con¬ 
served energy S is plotted for the disc mass A4 = 0.5M (these curves are 
given in solid diamonds; bottom axis applies). The top couple of curves has 
been obtained for the Paczyhski-Wiita potential, while the bottom (faster 
growing) couple for the logarithmic potential. In both cases, v grows al¬ 
most monotonously with A4 as well as with 8, having a single “dip” which 
is associated with the phase when the accessible region reaches above the 
disc edge. 


sen in Figs.lSttsfl 


v(^S^ jVl, disc) 


f — Veff(0 = 7r/2)] 27rrdr 

J 2'Kr dr 


(35) 


and plot the dependence of the result on M and £ for the PW and 
In potential in Fig.[T5] (Integration is performed over the accessible 
region; in cases where the the latter was not closed in the direction 
toward the centre, we have taken the lowest reachable r to be 5M.) 
In the ranges 0 < M < 0.8M and 0.945 < £+l< 0.965 (of 
which part is shown in the figure), and for both central potentials, 
the growth of v with either A4 or f is very similar. The comparison 
of plots shown in Fig.[T3thus suggests the following interpretation: 
the phase-space structure stays roughly the same for a moderate 
disc-mass perturbation, with the growing disc mass mostly induc¬ 
ing a shift of the orbits to higher kinetic energies. This aspect is 
surely present in the relativistic case studied in papers I-III as well, 
but would require a more subtle argument. 


deformation of geometry in the ring’s vicinity that the real physical 
distance (proper distance) to the ring comes out finite from out¬ 
side (when the ring is approached from bigger radii), whereas infi¬ 
nite from inside (when the ring is approached from smaller radii). 
When free motion is plotted in coordinates, the particles thus ap¬ 
pear repelled/attracted by the ring in the directions from which 
the ring is physically nearby/far away, i.e., they seem to be re¬ 
pelled towards larger radii, whereas attracted from smaller radii. 
The effect is strongest in the equat orial region. We noticed it and 
interpreted in Semerak et al.l ( Il999l) . and later this was repeated by 
IP’Afonseca et al. 1 2005 1. 

Since the above feature is “felt” up to several tenths of ring 
mass in the Weyl or Schwarzschild coordinates (in geometrized 
units), it might be somehow reflected in orbital statistics. How¬ 
ever, the effect is much better seen in the meridional plane (than 
in the equatorial one): the coordinate tracks of free particles, when 
approaching the ring from any latitudinal direction, are driven to¬ 
wards its inner side and hit it just along the equatorial plane. Inspec¬ 
tion of the ring’s neighbourhood in Poincare plots does not seem to 
indicate stronger anisotropy in the relativistic case. One can only 
observe slight differences in evolution of the main regular region 
centered just above the ring: for small ring mass, it is central sym¬ 
metric in all three descriptions, but when the mass reaches several 
percent of M, it “elongates” along the u” = 0 direction and fi¬ 
nally two new islands establish on its opposite radial sides, created 
by orbits circling around (“through”) the ring. This process starts 
somewhat before M = 0.02M in the relativistic system as well as 
in the system using the In potential, while in the PW-potential case 
it starts only before M = 0.03M. The only qualitative difference 
between the relativistic and the pseudo-Newtonian systems is that 
in the latter case, for large ring masses (from M = 0.8M for the In 
potential, while from M = 0.9M for the PW potential) a new pair 
of regular regions appear, again symmetrically with respect to the 
principal island, but now both lie above the ring radius. See mainly 
the last plot (A4 = I.IM) of the In-potential Fig. M where these 
two islands already dominate the section. It would be interesting 
to check whether the lack of this regular couple in the relativistic 
system has connection with the ring’s outward repulsion. 

However, it should be noted that the Poincare-surface analysis 
is best suited for the demonstration of long-term effects in the mo¬ 
tion of eternally orbiting particles, whereas the above mentioned 
feature mainly affects trajectories soon to be captured by the ring. 
Thus, the Poincare section will typically bear one or two points 
from such trajectories and their dynamical behaviour will be hardly 
discernible for most part. The only effect one could hope to observe 
in the surfaces of section is a deformation of invariant structures - 
of which we find no persuasive evidence. 


5.2 Remark on the Bach-Weyl ring 

The Bach-Weyl ring is actually an interesting source. Its poten¬ 
tial i29l is everywhere attractive, namely its field intensity (minus 
gradient of the potential) points toward the ring from all local lat¬ 
itudinal directions. In the Newtonian picture it thus represents an 
“ordinary” ring source. In relativity the potential remains valid, but 
the metric involves two functions, the second being given by a line 
integral of the potential gradient. In the BW-ring case, both func¬ 
tions are given by elliptic integrals and, as expected, both diverge 
at the very ring. The two divergences however combine to such a 

The “average speed” is certainly an ambiguous concept. We choose here 
a definition which is simple and natural. 


5.3 Resonance and chaos in orhit shapes 

Poincare surfaces of section represent a basic tool for assessing the 
overall structure of the possible test-particle motion, but one should 
keep in mind that they are really just sections through phase space, 
flattening out most of the information about individual trajectories. 
When comparing different systems, like the relativistic one and its 
pseudo-Newtonian counter-parts we are interested in here, one nat¬ 
urally first checks the Poincare diagrams for analogies and vari¬ 
ances, but in fact any statements concerning the occurrence of cer¬ 
tain structures in Poincare sections has to be taken with caution, 
because a particular sequence of recorded points (e.g. equatorial 
transitions) does not in general unveil a trajectory uniquely. 

In order to get an idea of what trajectory shapes such struc- 
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Figure 17. A counterpart of Fig.[T^ showing the shapes within the r sin 0, r cos 0 plane of the trajectories whose equatorial transitions have been recorded 
there (in the r, v'^ axes). The orbits are plotted up to some 200 periods and are coloured to be easily identifiable in Fig. [H From top left to right and bottom, 
the profile starts from the central orbit of the 3-fold island and proceeds toward the centre of the Fig. jl6l surface of section. All the plots have exactly the 
same scale, though the coordinate ranges (indicated along the axes in units of M) are adjusted to capture the orbits effectively. Orbits from “more interesting” 
regions are purple, one chaotic orbit is green. 


tures may represent and to illustrate what the statements about the 
frequency ratios mean for the actual trajectories, we select one of 
the sections obtained within the series capturing the dependence 
on iMMl-disc mass, namely the M = 0.35M section of Fig. 
|6] (where the black hole was simulated by the logarithmic poten¬ 
tial). This case represents the weakest perturbation for which sep- 
aratrix chaos already appears near the 3-fold island; the diagram 


is repeated in Fig. [T^ with a selection of orbits plotted in colours. 
The motion in the </> direction is dynamically unimportant (bound 
by conserved integral tj in the axially symmetric case, so we sup¬ 
press this dimension and illustrate the orbital shapes within the 
Weyl (p, z) meridional plane. The results are grouped in Fig. 1171 
marked by the same colours as their equatorial sections in Fig. M 

There are two distinct structures in Fig. [m the 3-fold and the 
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Figure 16. The Poincare diagram A4 = 0.35M from Fig. revisited 
with the aim to illustrate what kind of orbits its main structures represent. 
About 7400 transitions for each orbit has been recorded. The orbits are 
shown in colour to ensure their easy identification against Fig. [m where 
the meridional-plane shapes of their 200 periods are plotted. 


5-fold island; the ratio of the radial to vertical frequencies is 2:3 
for the former and 4:5 for the latter. The shapes of the trajectories 
reveal less thick resonances hidden in both the central and 3-fold is¬ 
land, but for most of them only a longer evolution track could con¬ 
firm whether it is actually a resonance or a near-periodic orbit only. 
However, one can notice a certain deformation due to the proximity 
of a resonance in Fig. [T3 the fish-like shape of the 2:3 trajectory 
corresponding to the 3-fold island is reflected in a significant part 
of the neighbouring non-resonant phase space, which might per¬ 
haps lead to observable signs in an ensemble of particles orbiting 
near the black hole. Besides obvious structures, one also notices, 
when recording data for Fig. d that the computation of the single 
purple orbit lying within the (blue) regular single-periodic region 
takes much longer time than that of the other orbits around. This 
typically indicates that one is close to a resonance, which is con¬ 
firmed in Fig. [it] Let us also point out to the rightmost orbit in 
the last-but-one row and to the middle one in the last row (both 
are purple): they are very similar, both lying just between a “box” 
regime and a resonant regime of the regular region, and analogous 
to a rotation of a pendulum very close to libration; the spatial cor¬ 
ridors more densely filled by the orbit in Fig.[T^then correspond to 
the pendulum near-stopping at the unstable top equilibrium (before 
falling back to rotation) which stands for an unstable counter-part 
of the stable periodic orbit at the resonance core. 

As seen in the second row of Fig. [T3 the time span corre¬ 
sponding to some 200 equatorial-plane intersections is not enough 
to discern between the regular trajectory (red) and the very close 
separatrix chaos (green). On the other hand, the respective surface 
of section in Fig. [T^ allows to discern order and chaos unambigu¬ 
ously at the toll of 7400 equatorial intersections. To better under¬ 
stand the computational/observational times required for a clear 
distinction between the regular and weakly chaotic orbit, we em¬ 
ploy a time-series recurrence analysis in the following section. 


5.4 Recurrence analysis 

ft is appropriate to support the Poincare-section observations by 
some other independent method. Like in paper 11, we turn to two 
recurrence methods here, one based on statistics over directions in 
which the orbits traverse a pre-selected mesh of phase-space cells. 


the other built on recurrences themselves to the neighbourhoods of 
phas e-space points. _ 

iKaplan & Glass! ( Il992h suggested to monitor the evolution of 
a tangent to the trajectory in small subsets of phase space which 
are crossed recurrently. For this purpose, the phase space is “recon¬ 
structed” from a given data series x(t) (either computed or mea¬ 
sured) by adding the latter’s replicas delayed by some shift At and 
its multiples. The method was designed to distinguish between de¬ 
terministic and random systems, but we saw in paper II that it is 
quite sensitive to weak irregularities and thus very well able to also 
recognize how chaotic the (deterministic) system is. Without going 
into details (see paper II for description of how we use the method 
for our system), let us only recall main points: 


• First the dimension d is chosen of the phase space to be re¬ 
constructed, plus the delay At and the size of boxes into which the 
phase space is divided. 

• Average tangents of a trajectory within a given (jth) box are 
summed (vector addition) for a large number of recurrent transits 
and the length of the result is suitably normalized; the result is de¬ 
noted as Vj (At). 

• The resulting norm is averaged then over all boxes which were 
crossed exactly n-times. 

• The result depends on n, on d, on At and on the lattice-box 
size. (The choice of these parameters in turn depends on how long 
data series one deals with.) With n it decreases roughly as 

for random data, whereas more slowly for a deterministic system 
(in a theoretical limit of an infinitely long series and infinitesimally 
fine grain, it even remains 1 for the deterministic case). The depen¬ 
dence on At is specifically studied on the deviation of the result 
from the value obtained for random walk, computed for each box 
and then averaged over all occupied boxes, 

' [Vj{At)]^ - (Rt- 


A = A(At) := 


i-(Rt)^ 


{Rn- is the average displacement per step for random walk of 
length rij in d dimensions). In a theoretical limit, A = 0 for a 
random walk, whereas A = 1 for a deterministic system; in prac¬ 
tice, A falls off roughly as autocorrelation function for a random 
series, while more slowly for a deterministic one. 


We have subjected to the Kaplan-Glass test two orbits from 
Figs.[T2and[T7] namely the outmost of the red-colour (3-fold) reg¬ 
ular ones and the nearby green-colour chaotic one which has arised 
from a separatrix breakup. The autocorrelation corresponding to 
the dependence of the “directional-vectors average” A on time de¬ 
lay At clearly confirms the different character of the orbits. Let us 
specify that we started the analysis from reconstructing the phase 
space as three-dimensional and dividing it in 25^ = 15625 boxes; 
average number of transitions through one box (among those which 
were crossed at least once) has been around 50. 

The second method rests on the statistics of recurrences to pre¬ 
scribed neighbourhoods of phase-space poin ts (either of the “orig¬ 
inal” phase space, or the reconstructed one). iMarwan et alj ( l2007h 
elaborated various useful outcomes of such a statistics and codified 
their computation in the RECURRENCE PLOTS software; we already 
applied it, in paper If, to the exact relativistic system. 

The main object of the analysis is the symmetric recurrence 
matrix 


Ri,3{e)=Q{e-\\Xi-Xj\\), i,j = l,...,N, (36) 

where Xi — X(Ti) denote N successive points of a given phase 
trajectory, e is the radius of a chosen neighbourhood (called thresh- 
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Figure 18. Two “neighbouring” orbits from Fig. [T6](and[T7j, namely the outmost of the red-colour (3-fold) regular ones (left plot) and the green-colour chaotic 
one (right plot), are clearly distinguished by the Kaplan-Glass “average directional vectors” recurrence method. The meaning of the A(At) dependence is 
explained in the main text. Recall that both orbits represent motion of free particles with £ + 1 = 0.955, £ = 3.75M in the field of a centre described by the 
logarithmic potential (mass M) and sun'ounded by the iMMl disc with mass M = 0.35M and inner radius rjisc = 20M. The orbits have been followed for 
about 500 OOOM of proper time (some 5000 periods); the top row shows the dependence A( At) from At = 0 up to At = 100 OOOM, while three selected 
intervals of At are added in more detail in the bottom row (the AT-axis labels are in thousands of M everywhere). 


old), 0 is the Heaviside step function and || ■ || denotes the chosen 
norm (we use a simple Euclidean norm, but other can be consid¬ 
ered, without significantly affecting the results). The matrix con¬ 
tains only units (meaning that j-th point is close to the i-th and rep¬ 
resented by black dots) and zeros (blank positions which mean the 
opposite). For regular systems, the recurrences arrange in distinct 
structures, in particular in long parallel diagonal lines and checker¬ 
board structures, whereas for random systems they are scattered 
without order; the chaotic systems provide something in between. 
A number of useful “quantifi ers” can then be extra cted from the re¬ 
currence data, as explained in iMarwan et al.l ( l2007h and also briefly 
reviewed in our paper II. The simplest ones follow from the lengths 
of the diagonal and vertical/horizontal lines which have occurred in 
the recurrence matrix. 

The recurrence pattern clearly depends on the time step At 
with which the trajectory is sampled and on the “target” radius e. 
Besides that, the matrix often contains false recurrence records that 
should be discarded from statistics. For example, if e is too large 
and the time step too small, several successive points of the trajec¬ 
tory of course lie within the e-neighbourhood of each other, but do 
not represent true recurrences. Due to the same reason, the real re¬ 
currences may then involve more than one point, even if the orbit 
comes across its certain previous part in quite a divergent manner. 
To overcome such false signals, several further “thresholds” are in¬ 


troduced and adjusted, mainly the minimal lengths of relevant di¬ 
agonal and vertical lines, (min and Umin- 

The choice of the recurrence threshold e should also take into 
account the physical extent of the orbit and its variance. However, 
we overcome this ambiguity by rescaling the time series in such a 
way that each variable has zero mean and unit variance. This also 
assures that motion in all coordinate directions have equal weight 
in the analysis, irrespective of the actual ranges spanned by the tra¬ 
jectory. 

For the relativistic-pseudo-Newtonian comparison using the 
recurrence-matrix analysis, we choose fig. 12 of paper II. There, 
several “quantifiers” were computed for 470 geodesics having 
specific energy S = 0.9532 and specific angular momentum 
I = 3.75M, sent tangentially (with u” = 0) from radii between 
r = 5M and r = 24M (with step 0.04M) from the equatorial 
plane of the system of a black hole (M) and the iMMl disc with 
M = 0.5M and raise = 18M. The orbits were followed for about 
250 OOOM of proper time with “sampling period” At = 45M, 
the minimal length of diagonal/vertical lines has been set at 90M 
and the radius of the recurrence neighbourhood (the threshold) at 
e = 1.25. Two of the quantifiers - the most simple one called DIV, 
given by reciprocal of the longest recurrence-matrix diagonal, and 
a much more “sophisticated” one read off from the slope of the 
histogram of diagonals (and providing an estimate of the maximal 
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Figure 19. Examples of the recun‘ence-plot results, obtained for free motion with £^(+1) = 0.9532 and i = 3.75M in the black-hole-disc field with 
Ad = 0.5M and r^isc = 18M. Exact relativistic system is represented in the left column, while pseudo-Newtonian system employing the logarithmic 
potential is in the right column. The top row shows Poincare diagrams coloured according to the value of DIV whose u'^ = O/u^ = 0 radial profile is also 
plotted in the middle row (going from blue to red, the value of DIV increases, which corresponds to increasing irregularity); the bottom row shows the same 
profile for another simple quantifier DET, given by ratio of the points which form a diagonal line longer than a certain minimum. The horizontal axes (r in 
units of M) are common for all rows and the vertical axes are common for both columns. One more remark: notice that in the left Poincare section the orange 
and red orbits are rather separated, whereas in the right one they are mixed within the chaotic sea. This is because the left section is actually divided into a 
“sticky” interior region harbouring weaker chaos and only slowly diffusing particles, and the outer chaotic sea with strong chaos. In the right section, the outer 
layer is mixed orange-red, with its less chaotic orange trajectories perhaps corresponding to motion “sticked”, for a short time, to the three small islands on the 
outskirts. For A4 ~ 0.65M the green layer of very weak chaos gets connected with the outskirts, thus yielding a picture rather similar to the relativistic case. 
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Lyapunov exponent), were particularly illustrated by colouring the 
computed orbits according to their values in the Poincare diagram. 
Two main observations were made: i) all the quantifiers proved sen¬ 
sitive to even tiny phase-space features, and ii) the computationally 
easy DIV quantifier proved equally efficient as the more sophisti¬ 
cated one. 

The above recurrence analysis was performed in a 6D phase 
space (r, 6, (j> and the respective velocities), while, for the present 
comparison, we have repeated it, for the same set of geodesics, 
in the (r sin6', r cos 6^) plane plus the respective velocity dimen¬ 
sions only. Elimination of (j) from the analysis has some interesting 
aspects, even though it is a Killing-symmetry coordinate. For in¬ 
stance, note that in the full 3D configuration space there are virtu¬ 
ally no true recurrences, since even the most regular central orbit is 
quasi-periodic in (f>. Within the meridional plane, on the other hand, 
the resonance cores produce true recurrences (see section [531 . 

Let us add that the (f) coordinate can be viewed as a kind of 
“dynamical memory”, because 



(37) 


(a relativistic formula only contains proper time r instead of t). 
Hence, the inclusion of (j) actually adds non-trivial information, so 
the change resulting from its elimination might indicate the robust¬ 
ness of various recurrence indicators. 

The parameters we have used for the re-analysis of fig. 12 of 
paper II are /min = 3 and Wmin = 3 for the minimal diagonal and 
vertical, the Theiler window w = 3 and the neighbourhood radius 
e = 0.8 (with Euclidean metric used for the distance). The trajec¬ 
tories for the recurrence matrix were then recorded at a time step 
At = 45M for a total time of about 250 OOOM like in paper II. 
As can be seen in Fig. the DIV indicator is not changed by 
the 3D—>-20 projection at all, whereas the DET quantifier turned 
out to be less robust in this respect. In particular, the DET quanti¬ 
fier seems to wrongly indicate that a large part of the central island 
is “less deterministic” than the surrounding chaos. To understand 
this point, let us recall that the DET indicator is defined as the 
ratio of the number of diagonal lines longer than /min to the num¬ 
ber of recurrence points. We checked that the orbits in the central 
island show a large number of recurrence points but not always 
grouped into longer diagonal lines. The performed normalization 
with respect to the total number of recurrence points thus has an 
undesirable effect in this case. 

Now to the comparison: we take an analogous pseudo- 
Newtonian situation, namely the gravitational system with “the 
same” parameters and with the central black hole simulated by the 
logarithmic potential (we do not employ the Paczynski-Wiita one, 
because that yields rather open accessible region, which makes the 
chaotic sea efficiently drained away to the centre), and subject it 
to the same recurrence-matrix analysis as performed in fig. 12 of 
paper II; the results are given in Fig. EH Clearly the phase-space 
structure is rich for the given parameters and also rather different 
from its paper-II counter-part. (As already stressed above, the rela¬ 
tivistic and pseudo-Newtonian systems are qualitatively similar, but 
the similar phase-space pictures are somewhat shifted with respect 
to each other in the parameter space.) The question is whether the 
corresponding recurrence patterns are still not alike, in spite of this 
first-sight difference. 

In Fig.[T3 the left column is relativistic and the right column is 
pseudo-Newtonian (with the logarithmic potential), both plotted in 
the same scale. Both Poincare sections are coloured by the longest- 
diagonal reciprocal DIV, the latter’s zero-velocity radial profile 


being also plotted below, and the last row shows another simple 
quantifier DET, given by ratio of the points which form a diago¬ 
nal line longer than a certain value within all the recurrence points. 
Although the surfaces of section reveal rather different structures, 
we do not see any big overall divergence in the recurrence charac¬ 
teristics. 


CONCLUDING REMARKS 

We have considered Newtonian dynamical systems describing 
the massive-test-particle motion in a gravitational field of a 
Schwarzschild-like centre simulated by a suitable potential and 
surrounded symmetrically by a gravitating thin disc or ring. Try¬ 
ing to learn how they differ from the corresponding relativistic 
system, namely the time-like geodesic dynamics in the field of a 
Schwarzschild black hole surrounded by “the same” disc or ring 
(described by the same Newtonian potential), we plotted Poincare 
diagrams of equatorial transitions for a number of similar situa¬ 
tions (same coordinate position and relative mass of the disc or 
ring, same values of the particles’ conserved energy and angular 
momentum) and found similar tendencies, typical for weakly non- 
integrable systems. The picture revealed by the surfaces of sec¬ 
tion was also confirmed by two recurrence methods, one resting on 
statistics over directions in which the orbits transit recurrently the 
boxes of a pre-selected phase-space mesh, and the other analysing 
the recurrences themselves to some prescribed neighbourhood of 
phase-space points. We have been using a different code than in 
previous papers of this series, so the present results also support 
robustness of the observations made. 

A careful conclusion would be that the pseudo-Newtonian ap¬ 
proach can reproduce the long-term dynamics of our relativistic 
system reasonably, though there appear various quantitative differ¬ 
ences. However, this conclusion strongly depends on the potential 
used to mimic the black-hole centre: the Paczynski-Wiita and the 
logarithmic potentials provide results very similar to the relativis¬ 
tic treatment, while the Nowak-Wagoner potential offers quite a 
different picture; some other potentials are not suitable for these 
purposes at all (although they may be efficient in another context). 
Yet even the Paczynski-Wiita and the logarithmic potentials dif¬ 
fer considerably (from the relativistic system as well as from each 
other) in the phase-space accessible region they determine - the 
PW potential is too open, whereas the In potential is too closed 
in direction toward the centre, which mainly affects how effec¬ 
tively the centre “sucks out” the chaotic orbits; nevertheless, this 
does not seem to influence much the behaviour of regular structures 
under parameter change. Generally, the pseudo-systems (involving 
the PW and mainly the In potential) can be labelled slightly more 
unstable than the exact relativistic system, since their phase-space 
structures evolve somewhat faster with parameters. 

As mentioned in preceding papers, there are several possi¬ 
ble directions of further study. One can certainly subject the dy¬ 
namical system - either the relativistic or the (pseudo-)Newtonian 
one - to still other methods (than already employed there), e.g. 
the Melnikov-integral calculation or the basin-boundary analysis, 
or to a more detailed study of its particularly “interesting” orbits 
(mainly the periodic ones). However, most important astrophysi- 
cally is to make our setting more realistic and to try to confront it 
with what is going on in real celestial systems. The simplest issue, 
at least within the static and axially symmetric situation, would be 
to add another gravitating components like central star cluster, a 
jet or a halo. Second, we plan to consider non-singular (i.e. 3D) 
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sources instead of the infinitesimally thin ones. This is especially 
important in systems where the orbits can reach very close to the 
sources, mainly in the relativistic description when the rings as well 
as edges of the thin discs usually represent a curvature singularity 
and (thus) the space-time is unnaturally deformed in their vicinity. 
In particular, one would be interested in how reasonably - as far as 
the long-term dynamics is concerned - the singular (Bach-Weyl) 
ring can approximate a toroidal source; this could be studied on a 
sequence of toroids gradually thinning to a ring. 

Once obtaining a sufficiently realistic description, it is desir¬ 
able to look for consequences for observational phenomena. For 
instance, the ensembles of initial conditions studied in the surfaces 
of section can be understood as actual collisionless ensembles of 
particles orbiting a black hole. The increased “suck-in” of the en¬ 
semble under perturbation then imply an enhanced accretion rate, 
while the resonances, on the other hand, correspond to particularly 
behaving oscillation modes. 

And then there are difficult aspects of “realisticity”: incorpo¬ 
rating {adequately) rotation of the gravitating bodies (which brings 
dragging effects into the relativistic systems) and possibly also back 
reaction resulting from the non-test character of the orbiter. 
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